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Black-hole-neutron-star mergers resulting in the disruption of the neutron star and the formation of an accre- 
tion disk and/or the ejection of unbound material are prime candidates for the joint detection of gravitational- 
wave and electromagnetic signals when the next generation of gravitational-wave detectors comes online. How- 
ever, the disruption of the neutron star and the properties of the post-merger remnant are very sensitive to the 
parameters of the binary (mass ratio, black hole spin, neutron star radius). In this paper, we study the impact 
of the radius of the neutron star and the alignment of the black hole spin on black-hole-neutron-star mergers 
within the range of mass ratio currently deemed most likely for field binaries (A/bh ~ 7Mns) and for black 
hole spins large enough for the neutron star to disrupt ( Jbh/A/bh = 0.9). We find that: (i) In this regime, 
the merger is particularly sensitive to the radius of the neutron star, with remnant masses varying from 0.3A/ns 
to O.IMns for changes of only 2 km in the NS radius; (ii) O.OIM0 — O.O5M0 of unbound material can be 
ejected with kinetic energy J> 10 51 ergs, a significant increase compared to low mass ratio, low spin binaries. 
This ejecta could power detectable post-merger optical and radio afterglows, (iii) Only a small fraction of the 
Advanced LIGO events in this parameter range have gravitational-wave signals which could offer constraints 
on the equation of state of the neutron star (at best ~ 1% — 2% of the events for a single detector at design 
sensitivity), (iv) A misaligned black hole spin works against disk formation, with less neutron star material 
remaining outside of the black hole after merger, and a larger fraction of that material remaining in the tidal tail 
instead of the forming accretion disk, (v) Large kicks ^kick <; 300 km/s can be given to the final black hole as 
a result of a precessing BHNS merger, when the disruption of the neutron star occurs just outside or within the 
innermost stable spherical orbit. 

PACS numbers: 04.25.dg, 04.30.-w, 04.40.Dg, 47.75.+f 



I. INTRODUCTION 



The next generation of ground-based gravitational-wave 
detectors (Advanced LIGO, Advanced Virgo and KAGRA Q|- 
Ufl) will progressively begin taking data over the next decade, 
opening an entirely new way to observe the universe. One 
of the sources of gravitational waves that they will de- 
tect are compact binary coalescences: binary black holes 
(BBH), black-hole-neutron-star (BHNS) and binary neutron 
star (BNS) systems J3]. In the presence of a neutron star, these 
gravitational-wave signals could be accompanied by electro- 
magnetic emissions, which would provide better sky local- 
ization and additional information about the characteristics 
and the environment of the binary. The most energetic, and 
most often discussed potential counterparts are short gamma- 
ray bursts (SGRBs, see e.g. JUl), while other possibilities in- 
clude the x-ray and optical afterglows of SGRBs, optical tran- 
sients due to the radioactive decay of neutron-rich unbound 
material, and radio emission from that ejecta as it decelerates 
in the interstellar medium (see Jit] for more details on EM 
signals emitted by compact binary mergers, and J7I for their 
detectability). The ejection of a small amount of material at 
ultrarelativistic speeds as a result of a shock in the region in 
which two neutron stars first get into contact was also recently 
proposed as a potential outcome of BNS mergers H. Finally, 
pre-merger electromagnetic transients are also a possibility, 



and could for example be due to the breaking of the neutron 
star crust |0|. 

Which of these effects occur in practice depends strongly 
on the parameters of the binary. The exact conditions leading 
to the emission of a SGRB are not known, and will depend 
on the physical process responsible for these bursts — most 
likely either the extraction of black hole rotational energy by 
the magnetic field of the accretion disk (Blandford-Znajek ef- 
fect BlOlO . magnetically-driven outflows in the accretion disk 
(Blandford-Payne effect Ullo . or the production of ultrarel- 
ativistic e + e~ pairs from the annihilation of vv pairs, them- 
selves produced by a hot accretion disk surrounding the rem- 
nant black hole J5I [HI]- In all cases, the presence of both 
an accretion disk and a baryon-poor region in which a rela- 
tivistic jet can be produced appears to be a prerequisite. The 
amount of matter required in the disk depends on the effi- 
ciency of the jet production mechanism and the energy of the 
burst, with estimates s panning multiple orders of magnitudes 
(lO _4 Af - O.5M ) lH Hi. The formation of an accretion 
disk is a natural result of the merger of BNS systems. For 
BHNS binaries, however, many initial conditions lead to the 
direct plunge of the neutron star into the black hole, before the 
tidal field of the hole can disrupt the neutron star and cause 
the formation of an accretion disk. Relatively massive black 
holes (A/bh ^ 7Mp)) are favored both by current popula- 
tion synthesis models 11151 [Trill and by the distribution of black 
hole masses measured in galactic X-ray binaries lfl7ll . In that 



2 



regime, tidal disruption only occurs for the most rapidly spin- 
ning black holes 1181 . 

The ejection of enough unbound material to power de- 
tectable electromagnetic transients is not a certainty either. 
General relativistic simulations of BNS systems have found 
that only a small amount of mass is unbound by the merger 
(M ej = 10" 4 M o - 1CT 2 M ) OH], while Newtonian 
smoothed particle hydrodynamics (SPH) simulations have 
more optimistic predictions (M e] ^ 10~ 2 M o ) JM]. BHNS 
systems have more asymmetric mass ratios, and are thus gen- 
erally more favorable for the ejection of neutron-rich mate- 
rial during the disruption of the neutron star. However, this 
requires the neutron star to be disrupted in the first place. 
Numerical simulations have shown that this will only be the 
case if the mass of the black hole is lower than what pop- 
ulation synthesis models predict, or if the spin of the black 
hole is high. The amount of material ejected in such mergers 
remains very uncertain. SPH results predict ejected masses 
of up to ~ O.IMq I2H1 . General relativistic simulations for 
low mass, low spin black holes find little ejected material 
(< O.OlM©) 12211 . but estimates for rapidly spinning black 
holes have not been offered yet. We show in this paper that 
more massive outflows (O.OIMq — 0.05Af Q ) are likely for 
high black hole spins. Finally, matter could be ejected due to 
magnetically-driven [23f] or neutrino-driven I124II winds in the 
disk. 

An important limitation of existing general relativistic sim- 
ulations of BHNS mergers is the lack of coverage of the range 
of black hole masses deemed most likely astrophysically. The 
only simulations considering mass ratios q = A/bh/A/ns > 5 
have shown that, even for very large neutron stars (i?NS ~ 
14.4 km) and aligned black hole spins, dimensionless black 
hole spins xbh = J/Mq^ <; 0.7 are required for tidal dis- 
ruption to occur lfl8ll . Simulations for more symmetric mass 
ratios have however already shown that smaller neutron star 
radii 1I22I l25ll or misaligned black hole spins f26il are likely to 
make tidal disruption harder. In this paper, we begin a more 
quantitative exploration of these effects. We consider BHNS 
binaries at mass ratio q = 7 and with large black hole spins 
Xbh = 0.9, and vary the radius of the neutron star and the 
orientation of the black hole spin. We currently limit our- 
selves to the most basic physical effects which will influence 
the dynamics of the merger: we treat gravity in a fully gen- 
eral relativistic framework, but use a simple T-law equation 
of state to describe the neutron star matter, and ignore the ef- 
fects of magnetic fields or neutrino radiation. More realistic 
equations of state are likely to influence tidal disruption, but 
to first order the compactness of the neutron star is expected 
to capture the main physical effects during merger |j22]. Mag- 
netic fields are important for the evolution of any post-merger 
accretion disk, as these disks are unstable to the magnetorota- 
tional instability (MRI), but very large internal fields are nec- 
essary for magnetic effects to influence the disruption of the 
neutron star Iu7l l28ll . Finally, while neutrinos are the main 
source of cooling in the disk, and are thus crucial to their evo- 
lution over a cooling timescale r„ (t v ~ 0.1s for the relatively 
dense disks considered by Lee et al. fT3l . and could be signif- 
icantly smaller for the lower density disks observed at the end 



of our simulations), neutrino emission will be negligible be- 
fore merger (neutron stars are expected to be extremely cool, 
with T « 1 MeV). This was confirmed numerically in the 
case of BNS 1I29I1 . Our simulations will thus capture the phys- 
ical effects which are important for the evolution of BHNS 
systems in the last few orbits before their merger, and during 
the merger itself. They are however not suitable for the long 
term evolution of the post-merger remnant. 

Gravitational waveforms from BHNS and BNS mergers are 
of particular interest for the constraints that they might offer 
on the unknown equation of state of the neutron star. Numeri- 
cal simulations indicate that from the last few orbits of a BNS 
merger occurring at 100 Mpc, constraints of ~ 1 km could be 
obtained on the radius of the neutron star ll30ll . Damour et 
al. Olfl have shown using Effective One Body waveforms that 
equations of state effects in the late inspiral could be measured 
for events of moderate signal-to-noise ratio (p ~ 16). Accu- 
rate numerical simulations are however necessary to calibrate 
such models at high frequency. Numerical results by Bernuzzi 
et al. l32tl have confirmed the predictions of Damour et al. ll3l"ll 
regarding the detectability of these equations of state effects, 
but existing simulations are not accurate enough to model the 
waveform with the accuracy required to take full advantage 
of all of the information that will be available in waveforms 
detected by Advanced LIGO. For BHNS mergers of nonspin- 
ning black holes, tidal effects during the inspiral are too small 
to be detected directly by Advanced LIGO iH HI]. But the 
cutoff in the gravitational-wave spectrum occurring when the 
neutron star is disrupted by the tidal field of the black hole can 
be. Semi-analytical models have been developed to attempt to 
extract that information ll34[|35ll . Numerical simulations map- 
ping the cutoff frequency across the relevant parameter space 
are however necessary to better calibrate them. 

At low mass ratios (q = 2 — 3) and for nonspinning black 
holes, Lackey et al. D6ll have shown from numerical sim- 
ulations that the combined effects of the tidal interactions 
during the inspiral and of the high-frequency cutoff of the 
signal would allow Advanced LIGO to detect variations of 
~ 10% — 40% in the radius of the neutron star for a favorable 
event at ~ 100 Mpc. This is thus slightly inferior to what can 
be done for binary neutron star systems located at the same 
distance from the observer (and, within a fixed volume, we 
expect many more BNS mergers than BHNS mergers). At 
higher mass ratio, tidal effects during the inspiral further de- 
crease ll33ll . On the other hand, due to the higher total mass 
of the system, the amplitude of the signal will be larger. Tidal 
disruption will also occur at lower frequency, and thus in a 
more favorable region of the LIGO noise curve. The spin of 
the black hole can also affect how much tidal distortion oc- 
curs during the inspiral (as the neutron star can get closer to 
a spinning BH before reaching its ISCO, tidal effects can be 
stronger). How these competing effects will affect our ability 
to measure the properties of neutron stars in BHNS binaries, 
or even to differentiate BHNS systems from BBH binaries, 
remains an important open question. 

In this paper, we study the influence of the radius of the neu- 
tron star and of the orientation of the black hole spin on the 
dynamics of the merger of BHNS binaries around the black 



FIG. 1: Numerical grid before disruption of the neutron star, below 
the equatorial plane of the binary. The black hole is the excised re- 
gion on the right (black sphere), while we superpose a linear color 
scale for the baryon density po- The subdomains on the outer edge 
of the plot connect to spherical shells covering the wave region. 



hole mass Mbh ~ 10A/© currently favored by population 
synthesis models, focusing on tidal effects during the inspiral, 
on the characteristics of the post-merger remnant (accretion 
disk formation, outflows), and on the properties of the emit- 
ted gravitational-wave signal — and in particular the effects of 
the neutron star radius on that signal, and the conditions under 
which such effects might be detected by the next generation 
of gravitational-wave experiments. We will begin by describ- 
ing briefly the numerical setup used for our simulations, as 
well as modifications to the code since the publication of Hill 
(Sec.|nJ. We will then detail the initial configurations evolved 
(Sec. ITTTb - and estimate the accuracy of the results (Sec. IIVb . 
Finally, the main physical results are presented in Sec. [V] 



II. NUMERICAL SETUP 

The numerical simulations presented here are performed 
with the SpEC code I37h . which evolves Einstein's equations 
of general relativity coupled to the relativistic hydrodynamics 
equations. Einstein's equations are solved using pseudospec- 
tral methods, in the generalized harmonics formulation 113811 . 
and excising the black hole interior. The numerical grid on 
which Einstein's equations are solved consists at first of 8 
spherical shells immediately around the black hole, 8 spher- 
ical shells and 1 inner ball in the region close to the neutron 
star, 24 spherical shells covering the far-field region, and a 
set of 13 distorted cylindrical shells and filled cylinders con- 
necting them (see Fig. [TJ. All subdomains are touching but 
not overlapping. The low, medium and high resolution runs 
correspond to a total number of points N = 57 3 , 64 3 , 72 3 . 

Once the neutron star disrupts, we cannot take advantage of 
an approximate spherical symmetry around the neutron star, 
and have to modify the pseudospectral grid. The shells around 
the black hole are replaced by a set of 264 distorted cubes, 
in order to allow for high angular resolution as the neutron 
star falls into the hole. The wave zone is still covered by 24 



FIG. 2: Same as Fig.fT] but during the disruption of the neutron star. 



spherical shells, while the region around the neutron star and 
the near field region are covered by non-overlapping distorted 
cubes (see Fig. |2). The resolution is chosen adaptively, by 
requiring that the relative truncation error for each set of ba- 
sis functions (measured from the coefficients of the spectral 
expansion of the evolved variables) is (0.5, 0.7, 1.0) x 10~ 4 
for our 3 resolutions. The actual number of grid points 
thus vary during the merger, peaking as the neutron star ac- 
cretes onto the black hole. At medium resolution, we have 

N ~ 70 3 - no 3 . 

The relativistic hydrodynamics equations are solved on a 
separate finite difference grid 1391 , using high-order shock 
capturing methods: we solve the approximate Riemann prob- 
lem on cell faces using the WEN05 reconstructor I140l442ll and 
HLL fluxes |43Tl . The grid covers only the region in which 
matter is present, and expands/contracts at discrete times as 
needed. Before disruption, we use N = 100 3 , 120 3 , 140 3 
points for the 3 different resolutions. During merger, we use 
instead 120 3 , 140 3 , 160 3 points. For configurations in which 
the black hole spin is aligned with the orbital angular momen- 
tum, we only evolve the region above the orbital plane and im- 
pose symmetry conditions on that plane (the number of grid 
points in the direction orthogonal to the orbital plane is then 
divided by 2). 

Both evolution schemes are at their most efficient when the 
system is quasi-stationary in the grid frame, and thus when the 
coordinates of that frame are chosen so that the grid is comov- 
ing with the binary. Additionally, to avoid having to impose 
unknown boundary conditions on the excision surface inside 
the black hole, we have to control that surface to ensure that 
all characteristic velocities of the generalized harmonic sys- 
tem are pointing out of the computational domain l44ll . Both 
objectives are attained through the use of a time-dependent 
map between the frame of the numerical grid and the iner- 
tial frame J45|. The map is the composition of a distortion 
of the region immediately around the black hole of the type 
r — > r + f(r) Ylim c im(t)Y lm (9, </>), which controls the size 
and shape of the excision surface (Y lm (8, <ft) are the spheri- 
cal harmonics functions), a translation keeping the black hole 
center in place, and a rotation and global scaling to follow the 
orbital evolution of the binary. More details on the methods 
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used to control the parameters of these maps in BHNS sim- 
ulations can be found in Foucart et al. (2012) ifTill . A more 
complete discussion of the rotation map for precessing sys- 
tems will be presented in Ossokine et al.(in preparation). 

Compared to the simulations performed in [fljll, we also 
fixed an error in the algorithm responsible for communicating 
source terms between the two numerical grids, which could 
cause the time stepping algorithm to be effectively of lower 
order than expected. The correction leads to reduced errors in 
the observed trajectories and in the phase of the gravitational 
waveforms (see Sec. UVBY During merger, we also modified 
the filtering of the variables evolved in the generalized har- 
monics system. For each dimension for which the spectral 
expansion uses a set of Chebyshev polynomials (i.e. any di- 
mension of a distorted cube, the radial direction for spherical 
and cylindrical shells and the vertical direction for cylinders), 
previous simulations filtered the time derivative of the evolved 
variables according to 



a« ltered = a l; 



(e-( 



1) 



(1) 



where a n is the n th coefficient of the spectral expansion of 
the time derivative on a given set of basis functions, and N 
is the total number of polynomials used in this set. Instead, 
we now simply filter the top mode of the expansion of the 
evolved variables at each time step. This much weaker filter- 
ing reduces the number of basis functions necessary to obtain 
a given precision, and allows the spectral adaptive mesh re- 
finement method to estimate the truncation error of the expan- 
sion more accurately (the number of basis functions in each 
dimension is chosen adaptively so that the truncation error of 
the spectral expansion is fixed, but that truncation error can 
only be measured accurately for unfiltered modes). 



III. INITIAL CONFIGURATIONS 

The initial data for these simulations is constructed as de- 
scribed in Foucart et al.(2008) ll4al . The constraints that Ein- 
stein's equations impose on the initial variables are solved 
in the extended conformal thin sandwich approximation pTil . 
under the assumption that the system is in equilibrium in a 
frame rotating at angular velocity ilo and contracting with ra- 
dial velocity — v r . The values of fio and v r determine the 
eccentricity of the system, and are chosen iteratively in order 
to minimize that eccentricity 114 811 . We go through two itera- 
tions of that procedure, starting from quasi-circular orbits (i.e. 
initial data with v r = and ilo chosen so that the initial mo- 
tion of the neutron star center has no radial component). The 
residual eccentricities at the end of the iterative procedure are 
e ~ 0.002 — 0.004. The free variables in the initial data (con- 
formal metric, extrinsic curvature) are the weighted superpo- 
sition of an isolated black hole in Kerr-Schild coordinates and 
of an isolated neutron star in isotropic coordinates, following 
the method developed by Lovelace et al.(2008) ll49ll for binary 
black holes. A more detailed description of the modifications 
required to apply this method to black-hole-neutron-star sys- 
tems is given in Foucart et al. (2008) j46ll . 



TABLE I: Initial configurations studied. All binaries have a mass 
ratio of 1:7, with the black hole dimensionless spin magnitude being 
^bh = 0.9. Obh is the angle between the rotation axis of the black 
hole and the initial orbital angular momentum of the binary, Cns = 
M-ns/Rns is the compactness of the neutron star, and R^ s 1,0 
the Schwarzschild radius of that neutron star assuming an ADM mass 
in isolation of IAMq. The orbital parameters are the eccentricity e, 
the initial orbital angular velocity times the total mass of the system, 
MQ(t — 0), and the number of gravitational-wave cycles before the 
peak of the wave amplitude, N cyc i cs (approximate numbers given for 
the precessing systems, in which mode mixing makes this variable 
ill-defined). 



Name 


e BH 


Cns 


M=lAM e 




N cycles 


e(t = 0) 


R12i0 


0° 


0.170 


12.2 km 


0.0413 


20.5 


0.004 


R13i0 


0° 


0.156 


13.3 km 


0.0413 


20.3 


0.003 


R14i0 


0° 


0.144 


14.4 km 


0.0413 


19.7 


0.003 


R14i20 


20° 


0.144 


14.4 km 


0.0412 


~ 18 


0.003 


R14i40 


40° 


0.144 


14.4 km 


0.0413 


~ 17 


0.004 


R14i60 


60° 


0.144 


14.4 km 


0.0415 


~ 14 


0.002 



Two series of initial configurations are considered in this 
paper (see Table |I), chosen in order to study separately the 
effects of the radius of the neutron star and of the orien- 
tation of the black hole spin on BHNS mergers at realistic 
mass ratios. All configurations consider a black hole of mass 
M BH = 7Af NS ~ 10M Q and spin xbh = ^bh/M| h = 0.9, 
where A/bh is the Christodolou mass of the black hole, J| H 
its angular momentum, and A/ns the ADM mass in isola- 
tion of a neutron star with the same baryon mass Af^ s as 
the neutron star in the binary. The neutron star is initially 
nonspinning. As our initial data does not exactly represent 
a BHNS binary in quasi-equilibrium, small transients are al- 
ways observed at the beginning of the simulations. After those 
transients, the mass and spin of the black hole are slightly 
modified (by ~ 0.01%). All simulations describe the neu- 
tron star matter as an ideal fluid with stress-energy tensor 
Tpv = {pa{]~ + € ) + P) u i-i u v + -P.9/jh an d use a T-law equation 
of state of index T = 2: 



P = 



P 
Po 



PoT 



(2) 
(3) 



where p is the baryon density of the neutron star material, P 
its pressure, e its internal energy, its 4-velocity, k a free 
constant and T a variable related to the temperature of the 
fluid (P t h = PoT is the thermal pressure in the fluid). To 
obtain the physical temperature T p hy S of the fluid from the 
variable T, we assume that the thermal pressure is the compo- 
sition of an ideal gas component and a black body component: 



T = 



+ f 



nT 4 



(4) 



2m„ po 

where m n is the nucleon mass, k the Boltzmann constant, and 
/ a function of T reflecting the fraction of relativistic particles 
in the gas (see |[50ll5lll for details). 

In the first group of simulations, we consider black hole 
spins aligned with the orbital angular momentum of the sys- 
tem, and modify the radius of the neutron star between R = 
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12.2 km and R = 14.4 km (for M NS = 1.4M ). This is done 
by modifying the value of the free parameter k in the equation 
of state of the fluid. We only consider this simple variation of 
the equation of state as we know that, to first order, the radius 
of the neutron star is the most important contribution to the 
dependance of BHNS mergers on the equation of state of the 
fluid ll22tl . The largest neutron star (case R14i0) is identical to 
the configuration studied in Foucart et al.(2012) 11811 . except 
that the initial separation is larger than in our previous work. 
As we will see, these parameters study most of the transition 
between BHNS mergers resulting in the formation of massive 
disks and those having nearly no material left out of the black 
hole a few milliseconds after merger. 

As we are considering T-law equations of state, it is worth 
mentioning that any of these simulations actually represent a 
continuum of systems, as they are invariant through the rescal- 
ing 

M' = K * M (5) 
R' = K*R (6) 
T' = K *T (7) 

where M is the mass scale of the binary, R the distance scale, 
T the time scale, and K an arbitrary positive constant. A 
more useful description of the initial conditions would thus 
use quantities which are also invariant under this rescaling, 
i.e. the mass ratio q = Mns/-MbHi the compactness of 
the neutron star Cns = Mns/Rns, and the dimensionless 
time t = T/M (in units in which G = c = 1). Realis- 
tic neutron stars of 1. 4Mp probably have radii in the range 
i?Ns ~ 9 km — 14 km f52||. with the most likely values being 
i?Ns ~ 11 km — 12 km 15311 . The values of Cns considered 
here are thus more likely to be found in neutron stars of ADM 
mass around or slightly below 1.4A/ Q , while probably unre- 
alistically low for very massive neutron stars (A/ns ~ 2Af Q ). 

The second set of simulations considers variations of the 
orientation of the black hole spin while maintaining the equa- 
tion of state fixed (using the larger neutron star with Cns = 
0.144). In terms of disk formation, this choice of equation of 
state is clearly optimistic, although not unrealistic, and thus 
provides an upper bound on the mass remaining outside of the 
black hole after merger. Moderate misalignments of the spin 
of the black hole with respect to the angular momentum of 
the binaries are an expected consequence of the kick that the 
supernova explosion is likely to impart to the forming neu- 
tron star. The actual distribution of the misalignment angle 
Obh between the spin of the black hole and the orbital angu- 
lar momentum of the binary is currently unknown, although 
@bh ^ 90° should probably be favored lfl5ll . Here we vary 
Obh between 0° and 60°, as such misalignments are both 
physically realistic and covering the range of parameters over 
which the properties of the final remnant vary significantly (at 
least for BHNS systems of mass ratio q = 7 with black hole 
spins xbh = 0.9). Misalignments are also often quoted as the 
angle ?7bh between the black hole spin and the total angular 
momentum of the system. For the systems considered here, 
6 B H = (20°, 40°, 60°), we have ?/ B h = (7°, 14°, 21°). 

All simulations begin at a coordinate separation d = 
7.44M, where M = M B h + Af NS is the total ADM mass 



of the system at infinite separation. This corresponds to an 
initial orbital angular velocity Mr2 or bit(£ = 0) ~ 0.041, 
or an initial gravitational-wave frequency /gw(* = 0) ~ 

235 Hz . Over the course of the simulation, the bi- 

\ Ai N s j 

nary will go through 7—10 orbits before merging. 



IV. ACCURACY 

The combination of spectral and finite difference methods 
used in our simulations can make it difficult to obtain strict 
error estimates: spectral methods are exponentially conver- 
gent in regions in which all variables are smooth, but only 
show polynomial convergence in the presence of discontinu- 
ities (such as at the surface of the neutron star or at a shock 
front). The finite difference methods used to evolve the equa- 
tions of relativistic hydrodynamics are second order in smooth 
regions, and first order at the location of shocks. As the region 
in which we get first order convergence should be of mea- 
sure zero, we expect at least second-order convergence as we 
increase the resolution of the finite difference grid. In prac- 
tice, we generally observe much faster convergence between 
the 3 resolutions considered here, particularly for quantities 
evolved on the spectral grid (e.g. trajectories, gravitational- 
wave signal,...). A conservative estimate of our error would 
thus be to assume second order convergence between our 
medium and high resolutions - the actual error being some- 
where between that value and the optimistic estimate obtained 
by simply looking at the difference between those two simula- 
tions. The ratio between these pessimistic and optimistic error 
estimates is ~ 3 for the simulations presented here. 



A. Final Remnant 

Of the characteristics of the final remnant listed in Table ITT1 
the parameters of the black hole (mass and spin) are the most 
accurate, with relative errors e 1 ^ <S 0.3% (i.e. differences of 
~ 0.1% between the medium and high resolutions). Global 
mass measurements (disk mass, tidal tail mass) are already 
less accurate, with e^Jass ^ 15% (O.OIMns difference mea- 
sured in the final remnant mass of the medium and high res- 
olution runs for configuration R13i0). Finally, the maximum 
density within the disk is only order-of-magnitude accurate: 
the distribution of matter within the disk remains fairly asym- 
metric and time-dependent at the end of the simulation, and 
variations of ~ 50% within ~ 1 ms should still be expected. 
Measurements of the mass of unbound material and the prop- 
erties of this ejecta have similar errors, and are discussed in 
more details in Sec. IV A 31 



B. Waveform Accuracy 

The phase accuracy of the gravitational waveforms in the 
non-precessing simulations presented here is about a factor of 
2 better than in our last set of simulations lfl8ll . even though 
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FIG. 3: Phase error in the dominant (2,2) mode of the gravitational- 
wave signal for simulation R13i0. We show the phase difference 
between our standard and high resolutions, both with (dashed red 
line) and without (solid black line) aligning the waves in phase and 
time. The dash-dotted green curve shows an estimate of the error in 
the extrapolation method (obtained by comparing extrapolation using 
different polynomial orders). 



the evolutions are ~ 2 — 3 orbits longer. This is most likely 
due to the correction of an error which effectively decreased 
the order of the time stepping method used in our simula- 
tions. Fig.[3]shows the phase difference between the medium 
and high resolution of the inspiral of simulation R13i0, both 
without any matching (i.e. by directly computing the phase 
difference between the output of the 2 resolutions), and after 
matching the waveforms over one period of the radial oscilla- 
tion of the orbit, choosing a time and phase shifts minimizing 
the difference between the two waveforms in the matching re- 
gion. The first method is the most direct assessment of the 
effect of numerical errors on the phase of the gravitational - 
wave signal, which are here of the order of a few tenths of a 
radian during the inspiral. The second method is more useful 
when comparing waveforms obtained in simulations starting 
from different initial conditions, and shows how different the 
waveforms would look to a gravitational- wave detector. Fig. [3] 
shows that, for matched waveforms, differences of the order 
of a few percents of a radian or less cannot be resolved by our 
numerical simulations. 

When considering waveform accuracy, numerical errors 
due to the discretization of the evolution equations are how- 
ever only part of the problem. Another potential source of 
error comes from extracting gravitational waves at finite radii, 
and then using polynomial extrapolation to obtain the wave- 
form at null infinity J54f|. An estimate of the error due to 
this process can be obtained by comparing the waveforms ob- 
tained using different polynomial orders for the extrapolation. 
Fig. [3] shows that this error is ~ 0.01 rad. Phase differences 
of the same order can also be due to the eccentricity of the 
binaries, at least for the eccentricities e ~ 0.002 — 0.004 con- 
sidered here. This can easily be seen from the oscillations in 
the phase difference between different configurations shown 
in Fig. QT| (the oscillations in the phase difference between 
simulations R13i0 and R14i0 are smaller that those between 



R12i0 and R14i0 because the radial oscillations of the first 
two cases happen to be nearly in phase at the beginning of the 
simulation). 

Adding all these sources of error, we can thus estimate that 
differences between numerical waveforms are large enough 
to be measured at our current accuracy only if, for matched 
waveforms, we have S(f> > 0.05 rad during inspiral. 

V. RESULTS 
A. Non-Precessing Binaries 

1. Inspiral : Tidal Effects 

Before the disruption of the neutron star, the main differ- 
ences between a BHNS inspiral and a BBH inspiral are due 
to the finite size of the neutron star, and its distortion un- 
der the influence of the tidal field of the black hole. The 
tidal distortion of the neutron star, and in particular its ef- 
fect on the gravitational-wave signal, has already been stud- 
ied in the Post-Newtonian framework. During the early in- 
spiral, Hinderer et al. ll33l found that for BNS systems the 
tidal effects would only be detectable by Advanced LIGO 
for the most favorable configurations (i.e. the largest neu- 
tron stars). Over the last few orbits, Damour et al. [E3ll find 
that tidal parameters would be detectable for BNS mergers 
of moderate signal-to-noise ratio (p ~ 16). But as these ef- 
fects are significantly smaller for more asymmetric mass ra- 
tios, the detection of tidal effects through gravitational waves 
is much more difficult for BHNS systems. A more detailed 
discussion of the detectability of the neutron star equation 
of state in our mergers is offered in Sec. IV A 5 1 - but from 
Fig. QT| alone, where we show the phase difference between 
our 3 non-precessing simulations, it is easy to see that up to 4 
gravitational-wave cycles before the peak of the gravitational- 
wave signal (/gw Si 500 Hz) the difference between these 
cases is not resolved numerically. This could however be due 
either to the fact that tidal effects on the waveform are ex- 
tremely small, or to a failure of the simulations to capture the 
tidal distortion of the neutron star properly. Accordingly, we 
need to test that the neutron star is tidally distorted during in- 
spiral, and that these tidal effects scale as expected. We com- 
pute the quadrupole moments of the neutron star 




(where p = y/gWpo, W = yl-j- g^UiUj, Ui is the spatial 
component of the 4-velocity and dV is the volume element, so 
that / pdV = M^g), and assume that they are due to first or- 
der to the composition of the tidal distortion of the neutron star 
and of a coordinate boost, acting along orthogonal directions. 
Qij is then diagonal in the coordinate frame (e±, e z ), where 
e± are two orthogonal unit vectors in the equatorial plane of 
the binary and e z is a unit vector in the direction of the orbital 
angular momentum (by symmetry, Q xz = Q yz = 0). The 
orientation of e± in the equatorial plane is a priori unknown, 
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and practically determined by solving for the rotation matrix 
diagonalizing . To first order, we assume that the tidal dis- 
tortion causes the neutron star to stretch along the direction 
e+ and contract along e_ and e z , while the boost causes a 
contraction along e_ and a stretch along e + and e z , i.e. 
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(10) 
(11) 



where Q is the magnitude of the tidal distortion, and B the 
amplitude of the boost distortion. From this decomposition, 
we can also retrieve the lag angle ^tidai between the tidal 
bulge and the line connecting the center of the hole and the 
center of the neutron star (i.e. the angle between e + and the 
line connecting the two centers). 

These quantities are clearly dependent on the coordinate 
system chosen. We cannot entirely remove that depen- 
dence, but can at least get a reasonable normalization for the 
quadrupole moments from the quantity 



It 



on 



pr 2 dV. 



(12) 



For all simulations, we find similar boost components 
B/Iqo ~ 2%: B/Iqo is a function of the binary separation, 
but does not depend on the equation of state considered. The 
lag angle is fairly constant too, with "ftidai ~ 20° — 25° 
for separations d ~ 60 km — 90 km ~ 4/A/bh — 6A/bh- 
The tidal component, on the other hand, varies strongly with 
both the binary separation and the equation of state. In 
Fig. |4] we show Q/Iqq for the three different equations of 
state considered here, and in the range of binary separations 
d ~ 60 km — 90 km. The tidal distortion goes from being of 
the same order as the boost effect at d ~ 90 km to a factor of 
2 — 3 larger at d ~ 60 km. Not surprisingly, the larger neutron 
star is significantly more distorted. 

To leading order, we expect the tidal distortion Q and the 
tidal field of the black hole (~ Mbh/c? 3 ) to be related by 



Q ~ 2k 2 R 
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(13) 



where k 2 is the tidal Love number of the neutron star (for 
r = 2 polytropes, which at T = are equivalent to the T- 
law equation of state used in our simulations, k 2 was com- 
puted by Hinderer HH). To verify that the tidal effects scale 
as expected, we thus compute the normalized tidal parameter 
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The superscript R14 refers to values for simulation R14i0, 
and do = 123 km is the initial separation of the binary. The 
normalization I^ 14 is the value of I 00 for simulation R14i0 
at the separation d at which we are measuring Q aoIm . The 



0.06 

^ 8 0.04 

0.02 


1.4 

1.2 

a 

1 

0.8 
0.6 



o R12i0 
o R13i0 
° R14i0 



60 



70 80 
d(km) 



90 



FIG. 4: Top: Tidal quadrupole Q normalized by Zoo- Bottom: Nor- 
malized quadrupole Qnorm (see eq. 114b . The dashed line shows the 
first order Post-Newtonian prediction Q no im = 1- 



numerical factor 0.0071 is computed in the limit d — > 00 
(i.e. with i^ 14 computed for an isolated neutron star), so 
that Q noTm (d — > 00) = 1. The scalings of k 2 , Cns an d d 
are chosen so that, as long as the tidal distortion of the neu- 
tron star follows the theoretical predictions, we will measure 

QnoTtn 1- 

In practice, our ability to measure Q norm accurately is lim- 
ited by the fact that the boost B and normalization Zoo only 
approximately model the distortion of the NS due to coordi- 
nate effects (i.e. the boost, but also any other gauge effect 
due to the coordinate choices made in the simulation). Mea- 
surements of Q„orm are thus unreliable for Q <, B ~ 0.02. 
Figure |4] shows that for Q ~ B, the scatter in the measure- 
ment of Qnorm is ~ 30%, while for Q ~ 2 — 373, it decreases 
below 10%. All measurements of Qnorm are compatible with 
Qnorm = 1 within that scatter, thus showing that the tidal 
distortion of the neutron star follows approximately the pre- 
dictions of Ref. 15511 . even at close separations. 

From these computations, we can thus conclude that the 
tidal distortion of the neutron star during the late inspiral is 
resolved by our numerical simulations, and scales with the 
binary separation and the equation of state of the neutron star 
in the manner expected from theoretical calculations, at least 
within SQ/Ioo - 0.005. 



2. Merger and Disk Formation 

An important question when considering BHNS mergers is 
the form of the post-merger remnant. To first order, this de- 
pends on whether the neutron star is disrupted before reaching 
the innermost stable circular orbit of the black hole or not. In 
the first case, a large amount of matter can remain outside 
of the hole after merger in the form of an accretion disk and 
a tidal tail. In the second case, no matter will remain. For 
a BHNS merger to be the progenitor of a short gamma-ray 
burst, the creation of an accretion disk is necessary. Accord- 
ingly, SGRBs are only possible if the neutron star disrupts. 
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TABLE II: Properties of the final remnant. M r 5 e ™ „ ant is the baryon 
mass remaining outside of the black hole 5ms after merger. M^ a is 
the baryon mass located at a coordinate radius greater than ~ 200 km 
at the same time. p„Iax is the maximum density in the disk, Xbh me 
dimensionless spin of the black hole at the end of the simulation, 
A/ BH the final Christodolou mass of the black hole and M the ADM 
mass of the system at infinite separation. 



Name 


Remnant 




p£S [io 11 


? /cm ] 


Xbh 






JU NK 


M 


R12i0 


0.10 


0.06 


2 




0.923 


0.960 


R13i0 


0.20 


0.11 


3 




0.919 


0.950 


R14i0 


0.30 


0.16 


21 




0.910 


0.935 


R14i20 


0.28 


0.15 


17 




0.909 


0.939 


R14i40 


0.15 


0.10 


3 




0.898 


0.959 


R14i60 


0.03 


0.03 


0.4 




0.862 


0.978 



Stellar disruption is facilitated by low black hole masses, high 
black hole spins and large neutron stars (see ||56ft for a simple 
fit to the results of previous numerical simulations). At low 
mass ratios (A/bh/^ns < 5), a moderately spinning black 
hole xbh ~ 0.5 is generally sufficient to provide disks of 
~ O.IMq. For the more massive black holes considered here, 
however, this is no longer the case. We have already shown 
that for spins xbh < 0.7 disk formation is unlikely even for 
large neutron stars (i?NS = 14-4 km). The simulations pre- 
sented here begin to explore how smaller neutron stars fare. 

Fig.|5]shows snapshot of simulations R14i0 and R12i0, the 
largest and smallest neutron stars considered here, both in the 
middle of the stellar disruption, and 5 ms later. The larger 
neutron star disrupts far enough from the black hole for a large 
portion of the matter to be initially ejected into a tidal tail, 
but the smaller neutron star disrupts just outside of the ISCO 
of the black hole. The top-right panel of Fig. [5] in particular 
shows how close the smaller neutron star is to the hole when it 
disrupts. From this picture, the fact than any material remains 
outside of the hole after merger is surprising in itself, and an 
indication of how strong the effects of the black hole spin can 
be on infalling material. 

Important differences are observed between the final state 
of these mergers. For the larger neutron star, 30% of the 
matter remains outside of the black hole 5 ms after merger. 
More importantly, about half that material has already formed 
a thick accretion disk, of about 100 km in radius and with 
peak density p^™* ~ 2 x 10 12 g/cm 3 . The formation of a 
hot, thick disk is less obvious for the smaller neutron stars. 
The amount of material remaining outside of the black hole 
is by no means negligible (10% — 20% of the neutron star, 
see Fig. |5J, but the maximum density is about an order of 
magnitude lower. In fact, if we look at the average surface 
density as a function of radius (Fig. |7), we see no evidence 
of an accumulation of higher density material at lower radii 
(~ 100 km), while that feature is clearly visible for the larger 
neutron star. From these results, we can also infer that smaller 
neutron stars i?Ns ~ H km would probably be unable to form 
any long-lived remnant. 

Evolutions including all the necessary microphysics (mag- 
netic fields, neutrino cooling) will be necessary to determine 
how these disks evolve over longer time scales (<; 0.1s). 



We can already see, however, that for all three configura- 
tions the material remaining outside of the black hole is hot 
(< T >~ 2MeV), and would be cooled by neutrino emis- 
sion. In that respect, the differences in density between the 
remnants could be significant, as they modify the opacity of 
the disk to neutrino radiation, and thus the efficiency with 
which the disk can transfer its energy into neutrinos. 

The properties of the final remnant are presented in Table HI] 
Comparing our results for the amount of material remaining 
outside of the black hole at late times with the predictions of 
Ref. HH], we find good agreement (within 2% — 3% of the 
neutron star mass) for the two smallest stars. The largest star 
forms a disk heavy enough that we are out of the range in 
which the predictions of Ref. |56|] are expected to be valid 
— and indeed, the disk formed in the simulation is signifi- 
cantly more massive that what Ref. ll56ll would predict. We 
also find consistency between our simulations and Ref. ll56ll 
on the neutron star radius below which no matter will remain 
outside of the black hole after merger (~ 10.5 km — 11 km). 
A more careful examination of the differences between our 
numerical results and Ref. j5oTl indicates that in the regime of 
high spin, high black hole masses considered here, the rem- 
nant mass probably has a steeper dependence on the radius of 
the neutron star than what would be guessed from Ref. ll56ll . 
a model fitted mostly to lower mass systems. However, these 
differences could also be explained by the expected variations 
in the remnant mass due to the internal structure of the neu- 
tron star (i.e. the fact that two neutron stars with the same 
radius but different internal structure will result in different 
post-merger disk masses), especially consideringthe fact that 
all of the simulations used to fit the model in I5a1 had larger 
tidal Love number &2 than the neutron stars from simulations 
R13i0 and R12i0. Overall, the magnitude of the differences 
between the numerical results and the model are roughly at 
the expected level. The final black hole masses and spins are 
also within the expected errors of existing analytical models, 
i.e. 1% — 2% away from the values derived by Pannarale ll57ll . 

For all configurations, the disruption and merger of the neu- 
tron star occurs over ~ 2 ms (see Fig. [6]). Mass accretion at 
later times is negligible compared with what is observed in 
lower mass ratio systems. When a disk forms, its main char- 
acteristics are however fairly similar to the lower mass ratio 
cases — except of course for the aforementioned lower den- 
sities and larger disk radii, which are a natural consequence 
of the higher black hole mass. Fig. Q shows a few character- 
istics of the forming accretion disk for the most strongly dis- 
rupted case. The surface density peaks at a distance of 100 km 
from the black hole, and the disk extends to about 150 km 
l . The orbital velocity profile is slightly sub-Keplerian (by 
about 10%), while the sound speed is ~ 0.25f2r and thus com- 
patible with a thick, thermally supported disk of scale height 
H = 0.25?*. This is consistent with the actual scale height of 
the disk, H ~ 0.2r — 0.3r. Finally, the inner edge of the disk 



Distances are measured in terms of the circumferential radius in the equa- 
torial plane, r c i rc = ^ f " y/g^dcp, where cj> is the azimuthal angle. 
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FIG. 5: Merger for the non precessing cases R14i0 (left) and R12i0 (right). The top panel shows the system at the time at which half of the 

x 10 10 g/cm 3 . The bottom 
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neutron star material has been accreted onto the black hole. We show densities down to p m i n 

panel shows the remnant 5ms later, plotting densities down to p m i n ~ 10 _8 Mq 2 ~ 6 x 10 9 g/cm 3 and cutting out the x > 0, y < quadrant. 
The differences in scale between the 4 figures can be determined knowing that the size of the black hole is always Rbh ~ 15 km. 
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FIG. 6: Baryon mass remaining outside of the black hole as a func- 
tion of time, for the 3 non-precessing cases R14i0, R13i0 and R12i0. 
We shift all the curves by the time t$o% at which half of the matter 
has been accreted onto the black hole. 



is particularly hot: we plot an estimate of the entropy 
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(15) 



(the effective constant K e g is defined by P = ft e ff/°o)> an( ^ ^ n< ^ 
s ~ 9 for r ~ 60 km. Within the disk, we still have s ~ 4— 5. 
As the disk settles down over ~ 10 ms — 20 ms, we would 



expect the entropy to exhibit a minimum at the peak of the 
surface density distribution, as was observed in lower mass 
ratio systems 0221 l26ll . 



3. Ejecta 

The ejection of unbound material by compact binary merg- 
ers is a prerequisite for some electromagnetic counterparts, 
most notably emissions due to the radioactive decay of the 
neutron-rich ejecta |@, [HI]. This ejecta can be obtained 
through various physical processes: unbound material in 
the tidal tail, but also magnetically-driven ll23ll or neutrino- 
driven j24ll winds. The study of winds goes beyond the scope 
of this article, as this requires accurate long-term evolution 
of the remnant disk and the inclusion of physical processes 
that are neglected in this work (magnetic fields, neutrino ra- 
diation). We will thus limit ourselves to the measurement of 
unbound material in the tidal tail. 

Even the presence of ejected material in the tidal tail can be 
difficult to assess in general relativistic simulations, particu- 
larly at high mass ratios. Maintaining a high enough resolu- 
tion in both the disk-forming region and the tidal tail is chal- 
lenging, and in practice matter can only be reliably evolved 
up to a distance of a few hundred kilometers from the black 
hole. This is indeed one of the main disadvantage of any grid- 
based simulations when compared with smoothed particle hy- 
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FIG. 7: Top: Average surface density 5ms after merger for the 3 non- 
precessing cases R14i0, R13i0 and R12i0. Bottom: Disk profile for 
simulation R14i0. Shown are the angular velocity (normalized by the 
velocity for circular orbits in this metric), the entropy and the sound 
speed (normalized by fir) as a function of the circumferential radius 
r. 



drodynamics methods, which can easily follow the evolution 
of tidal tails. In a time-independent spacetime and when pres- 
sure forces are negligible, it is easy to determine whether ma- 
terial is unbound: if Ut < — 1, then the material will escape to 
infinity (and — u t is the Lorentz factor of the fluid at infinity). 
This condition is also a fairly good approximation for low- 
density material far away from the central black hole after a 
compact binary merger, but becomes more and more inaccu- 
rate as one gets close to the black hole, or densities in the tidal 
tail become higher. 

One way to assess whether using the Ut < — 1 condition 
to find unbound material is accurate is to follow material over 
a sufficiently long period of time, and check that u t doesn't 
vary much. In our simulations, however, this only occurs for 
~ 1 ms before the material leaves the numerical grid, which 
leads to large uncertainties in the amount of unbound mate- 
rial, and its characteristics. A more detailed discussion of 
these issues will be presented in Deaton et al. (in preparation). 
Here, we limit ourselves to a discussion of measurements of 
Ut at relative low radii (<J 20Mbh ~ 300 km) and over short 
timescales, and note the uncertainties due to these approxima- 
tions. 

On Fig. [8] we plot the amount of mass with u t < —1 on 
the grid (and more than 60 km away from the black hole). 
The most compact neutron star, simulation R12i0, naturally 
has the least material in a tidal tail: about O.O9M , of which 
~ O.O15M appear unbound 2.5 ms — 3.5 ms after accretion 
onto the black hole begins (1.5 ms — 2.5 ms after the time 
at which 50% of the neutron star material has been accreted 
onto the black hole), with variations of only O.OO2M over the 
1 ms period over which measurements appear reliable. Mov- 
ing up in stellar radius, simulation R13i0 offers the most re- 
liable measurement of the ejected mass, with a stable value 
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FIG. 8: Mass of the unbound material, as measured by the condition 
u t < —1. For each configuration, t mcrgc is the time at which 50% 
of the neutron star material has been accreted onto the black hole. 
The dashed vertical line represent the time at which O.OOIM0 has 
escaped the grid (the low density tidal tail of simulation R12i0 cannot 
be followed accurately for more than 2.5 ms, at which point we stop 
measuring the mass of the ejecta). 



of M^- 13 ~ O.O46M long before matter starts flowing out 
of the grid (the total mass of the tidal tail, in this case, is 
~ O.16M ), and variations of 0.002A/ Q over 1 ms. From this 
run, we can also estimate the relative error due to finite nu- 
merical resolution in these measurements, and get e e j < 40% 
(at our "medium" resolution, we found M^ 13 ~ O.O4OM ). 
This is distinct from the uncertainties due to the approxima- 
tion made when using ut as a proxy for finding whether ma- 
terial is bound or not, and appears to be the dominant source 
of error for simulations R12i0 and R13i0. Finally, the largest 
neutron star shows the most uncertain measurements. Veloc- 
ities and densities in the ejecta are generally higher: the ap- 
proximate method takes more time to become accurate, but 
material remains on the grid for a shorter amount of time. 
There also is material with u t < —1 flowing directly out of 
the forming accretion disk, presumably as a result of shocks 
during disk formation, which makes it impossible to have all 
of the potential ejecta in the range 60 km < r < 300 km at 
any given time. Even by expanding the outer boundary of the 
grid by 50% compared with the two other runs, we thus find 
that the measured M^ 14 ~ 0.050Af Q only appear to settle 
at the time at which matter starts flowing out of the grid (and 
boundary effects might influence the properties of the ejecta). 
It is thus quite likely that the ejected mass is slightly larger 
than what is observed in the simulation. Overall, adding the 
two main sources of error (numerical resolution and use of 
Ut), we estimate that the ejected masses are 
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0.015M Q ± 0.010M o 
0.046M o ± 0.020M Q 
0.050M o ±0.035M Q . 



(16) 
(17) 
(18) 



From the measurements of u t , we can also determine the 
distribution of the velocity of the fluid at large distance from 
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FIG. 9: Distribution of asymptotic velocities for the unbound mate- 
rial of simulation R13i0. Two different times of the high-resolution 
simulation, and one time of the medium resolution simulation are 
shown. 
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FIG. 10: Gravitational-wave emission [(2,-2) mode] for the non- 
precessing simulations R12i0 and R14i0, extrapolated to infinity and 
normalized by the ratio D/M of the distance to the observer D to 
the total mass of the binary M. 



the black hole. This is shown in Fig. [9] for configuration 
R13i0. The distribution peaks at v/c ~ 0.2, and most of the 
ejecta has v/c < 0.5. The qualitative features of the veloc- 
ity distribution appear fairly robust when we vary the time at 
which we measure ut, and the resolution of the simulation — 
with more uncertainties for the high-velocity tail of the distri- 
bution. The merger with the more compact neutron stars has 
a very similar velocity profile. The situation is quite differ- 
ent for R1H0, where about half of the ejecta initially appear 
to have v/c > 0.5. At this point, however, we do not have 
the ability to follow such material for a long enough period of 
time to assess the reliability of the velocity estimates of that 
last configuration, and have to limit ourselves to the obser- 
vation that the ejecta appears more relativistic for the largest 
neutron star than in the other cases studied here. 

Finally, from these measurements, we can estimate the ki- 
netic energy of the ejecta, which would be available for fu- 
ture emission as it slows down in the interstellar medium. We 
should note that these results are only order of magnitude es- 
timates, as these energies are sensitive to the high-velocity tail 
of the velocity distribution, which is poorly constrained in our 
simulations. Additionally, our energy estimates are particu- 
larly unreliable for simulation R14i0, due to the large amount 
of poorly resolved high-velocity material that is rapidly leav- 
ing the grid. We find E e j ~ (1, 4, 40) x 10 51 ergs for simula- 
tions (R12iO,R13iO,R14iO) respectively. 

Even considering the large uncertainties in these measure- 
ments, it is interesting that we consistently find that in this 
region of the parameter space a few percents of a solar mass 
can be ejected from the system. This is indeed very different 
from the results obtained in the limit of low mass, low spin 
black holes (or for BNS lfl9lo . where only a negligible amount 
of material was found to be unbound. These results indicate 
that in the case of q ~ 7 BHNS binaries, tidal disruption of 
the neutron star (when it occurs) is likely to be accompanied 
by the ejection of ^ 10 _2 M Q of neutron-rich material. Such 
outflows are promising for optical emissions due to the ra- 
dioactive decay of neutron-rich elements in the ejecta, and the 



production of heavy elements resulting from r-process nucle- 
osynthesis. These ejecta might even be detectable as a radio 
afterglow as unbound material decelerates in the interstellar 
medium |@]: the kinetic energy available for radio emission is 
indeed larger than in supernova explosions. However, the lu- 
minosity of the radio afterglow heavily depends on the density 
of the environment, and BHNS mergers are likely to occur in 
much lower density environments than supernova explosions. 
The deceleration of the ejecta in the interstellar medium would 
then occur over longer timescales, and remain harder to de- 
tect. Massive ejecta in BHNS mergers are also limited to high 
spin configurations, and should not be considered as the norm 
unless xbh ~ 0.9 is standard for black holes in compact bina- 
ries. And the most energetic ejecta found here, in particular, 
is a very optimistic scenario for which, in addition of a high 
black hole spin, we used a very large neutron star. Finally, we 
should note that the detailed evolution of the tidal tail is likely 
to depend on details of the equation of state that our simple 
T-law model cannot capture. This problem should thus be re- 
visited carefully with a more realistic modeling of the neutron 
star fluid. 



4. Gravitational Waveforms 

Variations of the gravitational waveform emitted by BHNS 
mergers with the equation of state of the neutron star are 
mainly due to two effects: the small tidal distortion of the neu- 
tron star during the inspiral, which we discussed in Sec. IV A II 
and the cutoff in the gravitational-wave signal when the neu- 
tron star is disrupted (if tidal disruption occurs). In this sec- 
tion, we discuss measurements of these effects in our numer- 
ical simulations, while in the next section we will focus on 
their detectability by the Advanced LIGO detector. 

The effects of tides on the gravitational-wave signal of a 
BHNS binary before the disruption of the neutron star are ex- 
pected to be fairly small: Damour et al. ll3lll computed the 
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FIG. 11: Phase difference between the (2,-2) mode of the 
gravitational-wave emission of the non-precessing simulations. Both 
R12i0 and R13i0 are compared with R14i0. The dashed black lines 
show the matching interval, while the dashed green line shows the 
time at which the amplitude of the signal peaks for case R14i0 (the 
case in which the neutron star disrupts at the earliest time). 



phase difference 6^ due to tidal effects in the Fourier trans- 
form of the dominant [(2, 2)] mode of the waveform to 2.5PN 
order in the stationary phase approximation, 



S^ 2 T 5PN 



117 X 

8u 



x 5/2y2.5PN 



(19) 



where x = (Mwqw/2) 2 ' 3 is the standard Post-Newtonian 
parameter, wqw is the frequency of the (2, —2) mode of the 
gravitational strain, 



PN 



1 + 2.5a; 



-3/2 



8.51a; 2 - 3.927re 



5/2 



(20) 



and A is the tidal deformability parameter which for a BHNS 
binary is 



X 



l + 12g 



2fc 2 



26 3C* s (l + g) 



(21) 



For ugwM Si 0.2 and the binary parameters considered here, 
we get 6^ & 0.16 rad between simulations R12i0 and R14i0 
(or, for the phase <fi of the gravitational waveform in the time 
domain, 5(f) < 0.21 rad). However, these predictions have 
only been tested on BNS mergers up to x ~ Cns> while the 
tidal disruption of the neutron star in the binaries considered 
here occurs at X ~ 0.3 — 0.4. In that regime, the PN expansion 
no longer appears convergent: the 2.5PN term is of the same 
order as the leading order term. 

Measuring the small phase difference 5(f) in our simulation 
is a challenging problem. FigfTUl shows the dominant (2,-2) 
mode of the gravitational strain h(t) for simulations R12i0 
and R14i0, after the application of a time and a phase shift 
chosen to minimize the phase difference within a matching 
interval spanning 2 cycles of the radial oscillation frequency. 
Clearly, the waveforms are very similar up to the last cycle 
before the disruption of the larger neutron star. Fig. [TTlshows 
the phase differences between the waveforms of simulations 
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FIG. 12: Phase difference between the (2,-2) mode of the 
gravitational-wave emission of the non-precessing simulations as a 
function of the normalized gravitational-wave frequency lugwM. 
The dotted and dot-dashed curves correspond to the 1PN and 2.5PN 
predictions from Ref. j3lll . All curves are matched at ljgwM = 0.2. 
At lower frequency, residual eccentricity in the simulation makes 
measurements of 0(cjgw) too noisy to be useful. 



R12i0, R13i0 and R14i0. They remain under 0.05 rad until 
50 M before the peak of the gravitational-wave signal! More- 
over, the differences appear dominated by the influence of 
the residual eccentricity, not by equation of state effects. In 
Sec. IV A II we showed that the tidal distortion of the neu- 
tron star follows fairly closely Post-Newtonian predictions, 
yet this is clearly insufficient to have a measurable effect on 
the gravitational-wave signal for most of the evolution. 

A small phase difference between the waveforms can how- 
ever be hidden by the matching procedure used. Another way 
to compare the waveforms is to look at the phase <f> as a func- 
tion of the gravitational-wave frequency. This gets rid of the 
need to apply an arbitrary time shift to the waveforms. How- 
ever, the computation of wgw from the numerical waveform 
is fairly inaccurate, and even a small residual eccentricity can 
introduce a large noise in the resulting ^(wqw)- Computing 
the difference <5c/>(wgw) between two numerical simulations 
can thus only be done once the evolution of the orbital fre- 
quency due to orbital decay becomes fast enough to dominate 
the effects of eccentricity. Fig. [12] shows measurements of 
(5</>(cl>gw) for our numerical simulations. For ujq^M <, 0.2, 
there are large oscillations due to the eccentricity of the or- 
bit, and we cannot accurately measure <5</>(u;gw)- But in 
the frequency range 0.2 <, ujqwM <J 0.35, a phase shift 
of ~ 0.2 rad is clearly observed between simulations R12i0 
and R13i0, and the same between R13i0 and R14i0. This re- 
sult lies in between the predictions of the 1PN and 2.5PN ap- 
proximations. As the 2.5PN results are barely outside of the 
expected numerical error, these simulations are not accurate 
enough to improve on the PN predictions for tidal effects at 
high frequency. But we can confirm that the dephasing ob- 
tained from the 2.5PN predictions is at least correct within a 
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with / re f = 0.5 kHz (Note that any value of / re f in the 
range 0.3 kHz — 0.8 kHz gives nearly identical result as hf 7 / & 
is approximately constant during the inspiral, as shown in 
Fig. [T3T l. As the binding energy of these systems at t = 
is i?Q lnd = 0.0055M, we see that these binaries will radiate 
~ 2% — 2.5% of their energy before merging, and ~ 15% 
of their angular momentum (the more compact neutron stars 
naturally radiating more, as they disrupt later). The final kicks 
remain low (< 30 km/s — 50 km/s), as is generally observed 
for non-precessing BHNS binaries. 



J. DetectabUity of the Neutron Star Radius by Advanced LIGO 



FIG. 13: Spectrum of the gravitational-wave signal. h op t(f) is the 
spectrum of the dominant mode of the gravitational-wave signal as 
seen by an optimally oriented observer at 100 Mpc. The dashed 
curve shows the leading order PN behavior, h = Af~ 7 ^ , with 
amplitude A matched to the numerical results. The HighFreq noise 
curve of the Advanced LIGO detector is shown for tuning frequen- 
cies of 1kHz, 1.5 kHz and 2 kHz [See Eq. J26b: the noise curve 
tuned to 1 kHz is taken from Shoemaker et al. 115911 while the other 
two are approximate noise curves for different tuning frequencies]. 



TABLE III: Gravitational-wave emission over the course of the 
simulation as measured at a radius R = 275M, where M is the 
ADM mass of the system at infinite separation. E GW is the en- 
ergy contained in the waves, J GW their angular momentum, and 
«kick = F GW /M|H al the velocity kick given to that black hole. 
/St" is the cutoff frequency of the gravitational-wave signal defined 
byEq.[22] 



Name 


E uw /M 


J uw /M* 


Ukick(km/s) 


fcut (kHz) 


R12i0 


0.021 


0.16 


30 


2.1 


R13i0 


0.017 


0.15 


45 


1.8 


R14i0 


0.014 


0.13 


45 


1.5 


R14i20 


0.013 


0.13 


60 




R14i40 


0.013 


0.11 


150 




R14i60 


0.013 


0.09 


345 





factor of ~ 2 up to frequencies uigwM ~ 0.35 (/ ~ 1 kHz). 
The 1PN prediction also seem slightly more accurate than the 
2.5PN results. 

The effects of the disruption of the neutron star by the tidal 
field of the black hole are easier to see in the numerical re- 
sults. Fig. [13] shows the spectrum of the gravitational-wave 
signal for the three simulations R12i0, R13i0 and R14i0 (for 
an optimally oriented binary at 100 Mpc). The largest star 
disrupts earlier, and the gravitational-wave spectrum is cut 
slightly above 1 kHz. For the smallest star, the cutoff is at 
about 2 kHz. As opposed to tidal effects, this high-frequency 
cutoff is very well resolved in the simulations. 

Other physical quantities which can be extracted from the 
gravitational- wave signal are summarized in Table HH] the en- 
ergy and angular momentum radiated, the final kick given to 
the system, and a cutoff frequency / cut defined (arbitrarily) 
by the relation 



Keeping in mind the results of the previous section, we can 
begin to address another important question: the measura- 
bility of finite size effects on the gravitational waveform of 
BHNS mergers for mass ratios q ~ 7. An earlier analysis of 
these issues by Lackey et al. 13611 showed that at lower mass 
ratios (q = 2 — 3) and for nonspinning black holes Advanced 
LIGO would be sensitive to differences in the radius of the 
neutron star of order 10% — 40% for an optimally oriented 
BHNS merger located at 100 Mpc. This is due in part to the 
effects on the waveform of the tidal distortion of the neutron 
star, and in part to variations in the binary separation at which 
the neutron star disrupts and the gravitational-wave signal is 
cut off. 

At higher mass ratio, tidal effects are smaller. However, 
the disruption of the neutron star occurs at a lower frequency 
and the amplitude of the gravitational-wave signal is larger. 
It is thus unclear whether finite size effects will be easier or 
harder to detect. On Fig. [13] we show the spectrum of the 
gravitational-wave signal as seen by an optimally oriented ob- 
server located 100 Mpc away from the binary, and compare it 
with different Advanced LIGO detector's strain noise spectra 
(see below). At that distance the differences between the three 
simulations seem to be marginally measurable. In the rest of 
this section, we will attempt to quantify this statement more 
carefully. 

To determine whether the difference between two wave- 
forms h\ and h-2 can be detected by Advanced LIGO, we use 
the approximate condition llrjbll 



\\8h\\ 2 = (Sh, 6h) = (hi - h 2 , hi - hn) > 1. 
where the inner product is defined as 



(23) 



(9,h) 



' M ~9*{f)h{f) + ~g{f)h*{f) 
1 Sn(f) ■ ( } 



Here <?(/) and h(f) are the Fourier transforms of two wave- 
forms g(t) and h(t), and S n (f) is the one-sided power spec- 
tral density of the detector's strain noise, defined as 



/oo 
dre 2 " /T C n (r) , / > 0, 
-OO 



(25) 



2^(/cut)/I / t 6 



(22) 



where C„(r) is the noise correlation matrix for zero-mean, 
stationary noise. In this case, we will consider three of the 
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Advanced LIGO quideline noise curves defined in Ref. (5S 
the Zero Detuned Low Power spectrum, which is the expected 
sensitivity of the detector once signal recycling mirrors are in 
place, the Zero Detuned High Power spectrum, which is the 
final design sensitivity of Advanced LIGO, and a High Fre- 
quency noise curve optimized to take data at 1 kHz. We also 
approximate the effect of tuning the High Frequency noise 
curve to different frequencies by using the following estimate 
far £»(/): 



Sn(f) 



Q ZD-LP 



(/) s, 



(//«) 



ZD-LP 



(26) 



where S% D - LP (f) and S% F (f) are the Zero Detuned Low 
Power and High Frequency spectra from |p9il . a is a constant 
moving the tuning frequency to f — a kHz, and (3 is cho- 
sen so that / S n (f )df = J S^ F (.f)df. It is worth noting 
that these High-Frequency noise curves are for a detector op- 
erating with low power lasers, i.e. tuned versions of the Zero 
Detuned Low Power noise curves. 

When taking the inner product ||<5/i|| 2 , we choose one po- 
larization of hi and then allow for a time and phase shift in 
h,2, chosen to maximize the inner product {hi, hz) (which is 
identical to minimizing \\Sh\\ 2 for \\Sh\\ 2 <C ||^i,2|| 2 )- Fi- 
nally, we consider only the quadrupolar part of the waveform 
as measured by an optimally oriented observer: 



^opt 



(*) 



!(*)■ 



(27) 



In theory, the integration in Eq. (1241 should be carried 
over the entire frequency band of the detector. This is how- 
ever complicated for our waveforms, as the numerical simu- 
lations only cover the high-frequency part of the signal (/ > 
0.25 kHz). To construct a full waveform, the numerical re- 
sults should be hybridized with some analytical approxima- 
tion valid at low frequency (Post-Newtoninan, Effective One- 
Body,...). But in the region of parameter space considered 
here {q = 7, xbh = 0.9), the error coming from extend- 
ing these approximants to frequencies / ~ 0.25 kHz is sig- 
nificantly larger than the actual differences expected between 
waveforms for / < 0.25 kHz. Uncertainties in the construc- 
tion of the hybrid then dominate the measurement of ||<5/i|| 2 . 

We will instead consider three approximations to j| 6h\\ 2 : 

• Limiting the integration to frequencies / > 0.8 kHz, 
where the waveform is known exactly from the numer- 
ical simulations. This neglects tidal effects during the 
inspiral. 

• Limiting the integration to frequencies / < 0.8 kHz, 
and using PN predictions over the entire frequency 
range (see Appendix lAl. In this case, we ignore the ef- 
fects of the disruption of the neutron star, as well as er- 
rors in the PN predictions at high frequencies. As seen 
in the previous section, the 1PN predictions appear to 
remain fairly accurate for / < 0.8 kHz, and provide at 
least a good qualitative estimate of the tidal effects dur- 
ing the inspiral. 



TABLE IV: Detectability of tidal effects in non-precessing BHNS 
mergers for 3 different Advanced LIGO noise curves from Ref. l59ll : 
Zero-Detuned High Power, Zero-Detuned Low Power and High Fre- 
quency tuned at 1kHz (see also Fig.ll3t. The quantities \\Sh\\, de- 
fined in Eq. [23] are the detectability criteria for optimally oriented 
events at 100 Mpc. We consider first the results for the numer- 
ical waveform limited to / > 0.8 kHz (NR Only), then for the 
first order PN expansion [Eq. ([19}] limited to / < 0.8 kHz (1PN 
Only) and finally for hybrid waveforms matched in spectral space 
between 0.3 kHz < / < 0.8 kHz (Hybrid-IPN). For each case, 
we compare simulations [R12iO,R14iO] (||<5/i||i2-i4), [R12iO,R13iO] 
(||<Wi||i3-i3) and [R13iO,R14iO] (||<5A||i3-i4). For the hybrids, we 
also give results for approximate noise curves tuned to 1.5 kHz and 
2 kHz [see Eq. (26)]. 

NR Only 



S(f) 


\\Sh\\ 12-14 


\\Sh\ 12-13 


\\Sh |l3-14 


Zero Det H-P 


1.5 


0.7 


0.9 


Zero Det L-P 


0.7 


0.3 


0.4 


High Freq(l kHz) 


1.2 


0.4 


0.9 


1PN Only 


S(f) 


\\Sh || 12- 14 


||^||l2-13 


\\Sh 13-14 


Zero Det H-P 


1.0 


0.4 


0.6 


Zero Det L-P 


0.5 


0.2 


0.3 


High Freq(l kHz) 


0.4 


0.2 


0.3 



Hybrid-IPN 



S(f) 


\\Sh |l2-14 


||Jfc||ia-ia 


5^ 13-14 


Zero Det H-P 


2.0 


0.8 


1.2 


Zero Det L-P 


1.0 


0.4 


0.6 


High Freq(l kHz) 


1.4 


0.5 


1.0 


High Freq( 1.5 kHz) 


2.3 


1.0 


1.3 


High Freq(2.0 kHz) 


1.5 


0.8 


0.7 



• Combining the low-frequency PN predictions with the 
numerical results at high frequency by matching in the 
frequency domain the phase difference between two 
simulations to the predicted PN phase difference. This 
procedure is detailed in Appendix lAl 

The matching procedure and the differences between various 
PN orders each cause uncertainties of ~ 10% in \\Sh\\. 

The results for each method are summarized in Table [IV] 
We find that for the Zero Detuned noise curves, the high fre- 
quency cutoff is only ~ 50% more important than the low- 
frequency tidal effects. The tuned High Frequency noise 
curves, quite naturally, are much more sensitive to the dis- 
ruption of the neutron star than to tidal effects. For the Zero 
Detuned High Power noise curve, the detectability criteria is 
satisfied if the neutron stars have radii differing by 



D 



100 Mpc 



km, 



(28) 



where D is the distance to the observer. The Low Power noise 
curve requires the binary to be about twice as close. But re- 
sults slightly better than for the Zero Detuned High Power 
noise curve can be recovered by tuning the low-power detec- 
tor to observe at ~ 1.5 kHz. 

Differentiating a binary black hole system from a BHNS 
binary would only be slightly easier. The phase difference be- 
tween a point-particle waveform and the waveform of simula- 
tion R14i0 is, during inspiral, only ~ 1.5 times larger than the 
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phase difference between simulations R12i0 and R14i0. And 
the high-frequency cutoff in the waveform would not help us 
much, as the most compact star considered here (R12i0) does 
not disrupt very far from the ISCO. 

These results are not very promising. The condition 
\\5h\\ > 1 will only be satisfied for binaries with signal-to- 
noise ratio p > 30 — 40, or about 1% — 2% of the Advanced 
LIGO events (for the waveforms from simulations R12i0 and 
R14i0, and assuming that BHNS binaries are equally dis- 
tributed in volume). The criteria \\5h\\ > 1 is also an opti- 
mistic limit: it does not take into account the fact that all other 
parameters of the binary are here assumed to be known. Un- 
certainties in the masses of the objects and the spin of the 
black hole will significantly affect these results, making it 
harder to detect equation of state effects. Additionally, both 
the Post-Newtonian dephasing and the variations in the high- 
frequency cutoff of the signal are helped by the fact that we 
are considering a rapidly rotating black hole. For a nonspin- 
ning black hole, the neutron star would reach the ISCO at 
lower frequencies, and tidal effects during the inspiral would 
be even smaller. And as the neutron star would plunge into 
the black hole before being disrupted, there is no guarantee 
that there would be a measurable difference in the cutoff fre- 
quency of the waveform (although this question certainly de- 
serves further investigation: the merger would also occur in a 
more favorable frequency range). Finally, any real data anal- 
ysis would require the knowledge of the waveform at low fre- 
quency, which is not at this point known with enough accu- 
racy for binaries with q=-l and xbh = 0.9. The theoretical 
detectability conditions considered here are thus certainly too 
generous. 



B. Precessing Binaries 

Our second set of BHNS mergers considers variations of 
the orientation of the black hole spin. The starting configura- 
tion is the largest neutron star studied in the previous section. 
The black hole spin is inclined with respect to the orbital an- 
gular momentum by 6>bh = 20°, 40°, 60°. In all cases, the 
misaligned component of the spin is initially along the line 
connecting the black hole and neutron star centers. We do 
not expect this choice to affect the qualitative feature of the 
merger (disruption, disk formation) [261 . It should, however, 
influence the magnitude of the velocity kick given to the final 
black hole as a result of gravitational-wave emission lIcSlll . 



1. Inspired: Orbital Precession 

During inspiral, the main difference with the aligned con- 
figurations will naturally be the precession of both the black 
hole spin and the orbital angular momentum around the total 
angular momentum of the system. A simple coordinate mea- 
surement of that precession is presented in Fig. [14] using the 
angle 6*bi n between the direction of the initial orbital angu- 
lar momentum and the normal to the orbital plane (defined as 

(CBH - C N s) X («BH - VNS))- 
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FIG. 14: Inclination #bm of the line connecting the BH and NS center 
with respect to the initial orbital plane, up to the beginning of stellar 
disruption. 



We see that the amplitude of the precession of the orbital 
plane vary from 27° for R14i20 to 82° for the most inclined 
case. If the total angular momentum of the system was con- 
stant (i.e. in the absence of gravitational-wave emission), we 
would expect the amplitude of that precession to be twice the 
initial angle between the orbital angular momentum and the 
total angular momentum, i.e. 25° and 78° for simulations 
R14i20 and R14i60 respectively. Somewhat larger values are 
expected for radiating systems, as the loss of angular momen- 
tum due to gravitational-wave emission is to first order aligned 
with the current orbital angular momentum of the binary: tak- 
ing into account the angular momentum lost to gravitational 
waves listed in Table ITTT1 would correct these estimates to 30° 
and 85° (which are now overestimates, as part of the radi- 
ated angular momentum is emitted during the disruption of 
the neutron star). Over the course of the binary evolution, the 
most inclined binary goes through slightly more than half of 
a precession period, while R14i20 gets close to completing a 
full precession period (as #bin is defined with respect to the 
initial orbital plane, we get 6*bi n ~ after a full precession 
period, and not after half a precession period). 

A simple comparison of the measured precession of the or- 
bital plane with the Post-Newtonian predictions from rt62ll is 
also presented on Fig. [14] The IPN and 2PN curves are ob- 
tained by evolving the initial BH spin and orbital angular mo- 
mentum using the Post-Newtonian formulae from ll62ll . but 
assuming that the trajectories of the compact objects are those 
observed in the simulation (i.e. when computing the relative 
position and velocity of the objects, we use our results and 
not a Post-Newtonian evolution of the initial conditions, as 
the Post-Newtonian equations of motion are quite inaccurate 
so close to merger). We see that the period over which the 
orbital plane precesses in the simulations matches the PN pre- 
dictions very well, while the amplitude of the precession is 
~ 10% smaller in the numerical results. Such differences (as 
well as the additional oscillations in the value of the inclina- 
tion angle) could however easily be due to the fact that these 
measurements are certainly not gauge-independent. 
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A large precession of the orbital plane is quite natural for 
the high mass ratio, high spin systems considered here. In- 
deed, the angular momentum of the black hole is Jbh/M 2 = 
0.689 while the initial orbital angular momentum is only 
L or hit/M 2 = 0.354 (and about 10% of the total angular 
momentum, or ~ 30% of L rbit, is radiated in gravitational 
waves). This means that even for relatively low inclinations 
of the black hole spin (such as in simulation R14i20), large 
variations of the orientation of the orbital plane occur. These 
oscillations and, more importantly, the shift in the phase of the 
gravitational waveform with respect to aligned spin templates 
which accompany them, can make the detection of precessing 
binaries challenging. The higher dimensionality of the param- 
eter space to consider also complicates parameter estimates 
— although there are also positive effects due to precession: 
some of the degeneracies existing between the parameters of 
the system for aligned binaries are broken for precessing sys- 
tems l3ltl . As mentioned in the previous section, an inaccurate 
determination of the parameters of the system increases the er- 
ror in any measurement of the neutron star equation of state 
from gravitational waveforms. Obtaining proper constraints 
on the parameters of a BHNS binary with misaligned spin, 
which requires reliable templates for precessing systems, is 
thus a prerequisite to any attempt at constraining the equation 
of state of neutron stars from BHNS waveforms, at least for 
the high-spin configurations considered here (unless the spin 
of the black hole is aligned with the orbital angular momen- 
tum of the system by some unknown mechanism during the 
pre-merger evolution of the binary). 



2. Disruption and Disk Formation 

The effect of spin misalignment on the properties of the 
remnant of a BHNS merger were first studied in a general 
relativistic framework in Ref. ll26"ll . Material on an orbit in- 
clined with respect to the equatorial plane of the black hole 
reaches the region in which stable orbits no longer exist at 
a larger separation than material in the equatorial plane on 
a corotating orbit. Effectively, this means that BHNS merg- 
ers with misaligned black hole spins are roughly equivalent 
to mergers with a lower black hole spin aligned with the or- 
bital angular momentum. Disruption becomes less likely for 
misaligned configurations, and the mass remaining outside of 
the black hole at late times is smaller. The smallest radius 
at which stable orbits exist for black holes with xbh = 0.9 
at inclination 6> B h = (20°, 40°, 60°) is equal to the radius of 
the innermost stable circular orbit of a black hole with aligned 
spin xbh = (0.89,0.80,0.62) 00. 2 And indeed, the 
disruption of the neutron star is close to what is expected for 



2 The method used here to approximate the "effective" spin of a BHNS 
binary for the purpose of tidal disruption was first pointed out to us by 
Nicholas Stone. This method has been used to impose constraints on sys- 
tems which could lead to SGRBs in Stone et al. I63I . Our numerical sim- 
ulations tend to confirm that this is indeed a reasonable approximation to 
the result of tidal disruption in precessing BHNS binaries. 



such spins: for the low inclination case R14i20, the mass of 
material remaining outside of the black hole after merger is 
nearly identical to what was observed in the aligned config- 
uration R14i0 while in simulations R12I40 and R14i60, 15% 
and 3% of the neutron star mass remain outside of the hole 
5 ms after merger. This is very similar to what our simple fit- 
ting model Ipql would predict (13% and 1% of the neutron 
star matter surviving the merger for xbh = (0.80, 0.62)), 
or what might be inferred from numerical simulations for 
smaller black hole spins d (which found 6% of the mat- 
ter remaining outside of the hole for xbh = 0.7 and none for 
Xbh = 0.5, all other parameters being identical to the cases 
considered here). It thus seems fairly likely that, for the pur- 
pose of measuring the mass of material remaining outside of 
the black hole at late times at least, the disruption of the neu- 
tron star in BHNS mergers with misaligned black hole spins 
can be modeled with good accuracy by considering the results 
of aligned configurations only. 

However, there are important qualitative differences be- 
tween the behavior of aligned and misaligned configurations. 
Fig.[l5lshows snapshot of the two simulations with the lowest 
inclination angle for the black hole spin (R14i20 and R14i40), 
at a time at which 50% of the neutron star has been accreted 
onto the black hole (top) as well as 5 ms later (bottom). These 
can be directly compared with Fig. |5]for non-precessing sys- 
tems. Simulation R14i20 mostly behaves as the aligned case, 
although the slight precession of the tidal tail with respect to 
the disk induces small differences in the formation of the disk, 
and should affect its subsequent evolution as matter falls back 
from changing directions. At moderate inclinations (R14i40), 
the changes are more drastic. There is not much of a disk 
forming: most of the remaining material is in a long tidal 
tail, which is differentially precessing. For aligned configu- 
rations, a disk generally starts to form as the front edge of 
the accretion flow wraps around the black hole and hits the 
material from the tidal tail, creating a shock, outflows, and a 
rapid redistribution of the tidal tail material. Disk-tail interac- 
tions are much less visible in inclined simulations (although 
there are still contacts between the inner and outer edge of the 
tidal tail as it wraps around the black hole). Existing shocks 
are however still sufficient to heat the remaining material in 
simulations R14i20 and R14i40 to an average temperature 
< T 2McV - 3MeV. Not surprisingly the remnant 
of the most inclined merger (R14i60), which does not form a 
disk, is much cooler (T < 1 McV). 

The significant asymmetry of this system with respect to 
the equatorial plane of the black hole, as well as the differ- 
ential precession between fluid elements, could also have im- 
portant consequences for the long term evolution of the sys- 
tem. Etienne et al. j66ll have shown that asymmetries and mo- 
tion across the equatorial plane of the black hole are impor- 
tant in magnetized BHNS mergers as they contribute to the 
formation of a larger toroidal field within the disk (aligned 
configurations with equatorial symmetry, on the other hand, 
form disks with mostly poloidal magnetic fields, as the mag- 
netic field comes from the winding of the original field lines 
frozen within the disrupting neutron star). The toroidal field 
is amplified by the fastest growing mode of the magnetoro- 
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FIG. 15: Same as Fig. [5] but for the precessing binaries R14i20 (left) and R14i40 (right). 



tational instability(MRJ), the mechanism most likely to allow 
the magnetic field in the accretion disk to grow up to levels 
at which jets (and short gamma-ray bursts) could be created. 
Misaligned black hole spins could thus lead to qualitative dif- 
ference in the evolution of magnetized disks. They would also 
make it easier to resolve the growth of the MRI numerically, 
as the wavelength of the fastest growing MRI mode scales 
with the magnitude of the toroidal component of the magnetic 
field ll67ll . On the other hand, the disks formed from precess- 
ing BHNS binaries have lower densities and lower masses. 
Compared to aligned configurations, a larger fraction of the 
matter remaining outside of the black hole is sent in the tidal 
tail. The most inclined configuration studied here (R14i60), 
which barely disrupts, actually keeps less than 1% of the neu- 
tron star material within 200 km of the black hole. Nearly all 
of the remnant mass (3% of the neutron star material) is on 
highly eccentric, differentially precessing orbits. 

If the velocity kicks given to neutron stars during super- 
nova explosions are such that a majority of BHNS systems 
are within the range of inclination of the black hole spin stud- 
ied here (as might be expected from the results of Belczin- 
sky et al. Il5ll ). the effects of inclination on the orbital evo- 
lution of the binary, and consequently on the gravitational- 
wave signal, would be significant. But the conditions required 
for the disruption of the neutron star to occur in a significant 
number of systems would not be dramatically modified from 
those derived in non-precessing systems: the location of the 
marginally stable orbit does not change much for inclinations 



#bh < 30° (see e.g. HH), at least for the relatively large 
spins xbh > 0.7 which we already know are necessary for 
tidal disruption to occur for mass ratios q ~ 5 — 10. The 
remnant disks would be less massive, but with larger toroidal 
magnetic fields initially. Considering that even low mass disks 
(~ O.OIMq) are energetically sufficient to power gamma-ray 
bursts if the conditions (temperature, magnetic fields) are oth- 
erwise right 11131 Il4ll . misaligned configurations could in the 
end prove more favorable than the aligned cases. 



3. Waveforms 

The waveforms from precessing BHNS binaries are sig- 
nificantly different from their non-precessing counterparts, 
as previously mentioned: the precession of the orbital plane 
shown in Fig. [14] will cause a modulation of the preferred di- 
rection for gravitational-wave emission (see Fig. [16}, which 
has to be properly modeled in order to avoid significantly re- 
ducing our sensitivity to any waveform emitted by a precess- 
ing system. The modeling of precessing waveforms goes be- 
yond the scope of this article. A more detailed study of the 
impact of precession on the detectability of BHNS mergers 
can be found in Brown et al. 16 811 . Assuming an isotropic 
distribution of black hole spin, about half of the BHNS bi- 
naries within the theoretical range of the next generation of 
gravitational-wave detectors could be missed given the per- 
formance of the current search methods in the regime of high 
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FIG. 16: Gravitational waveforms as measured by observers who, at 
t = 0, see the binary face-on (solid black line) and edge-on along 
the line connecting the centers of the two objects (dashed red line), h 
is normalized by the ratio D/M of the distance to the observer D to 
the total mass of the binary M. We show both the h+ and h x polar- 
ization, with h x shifted by 0.2. Top: Non-precessing configuration 
R14i0. The face-on observer is always optimally located, while the 
edge-on observer receives a weaker signal in h+ (in which higher 
order modes are however more visible), and no signal at all in h x ■ 
Bottom: Precessing configuration R14i60. The envelope of the sig- 
nal varies in time as the binary precess, and the optimal orientation 
is modified accordingly. The merger also occurs at an earlier time, as 
the orbital hang-up is only due to the aligned component of the black 
hole spin. 



ied Ir26ll . This is due to the fact that most of the kick is received 
at the time of merger, while in BHNS systems at lower mass 
ratio the neutron star disrupts before merger, effectively cut- 
ting off the asymmetric gravitational-wave emission responsi- 
ble for the kick. Here, however, the disruption of the neutron 
star occurs very late in the inspiral — or nearly not at all in 
the case of the most inclined configuration. The strong dis- 
tortion of the neutron star as it plunges into the black hole 
will still reduce gravitational-wave emission at merger, even 
for configurations in which no matter remains outside of the 
black hole afterwards, but not nearly as much as for aligned 
configurations, or for more symmetric mass ratios. Accord- 
ingly, we find that much larger kicks i>kick ~ 345km/s are 
now possible. It should also be noted that the kick obtained 
from binary mergers is proportional to cos (4> + <fio), with (ft 
the orbital phase at merger and </>o an unknown phase shift. 
The only way to measure the maximum kick from a specific 
configuration is thus to consider a sequence of mergers, span- 
ning a range of phases </>. All that can be said from a sin- 
gle configuration is that kicks larger than 345 km/s are possi- 
ble. Without more studies of the dependence of «kick in <fi, we 
cannot in principle exclude the possibility that kicks for these 
configurations could be nearly as large as in binary black hole 
mergers, although maximal values closer to those presented 
here appear more likely. 



VI. CONCLUSIONS 



mass ratio, strongly precessing systems. Attempts to reduce 
the complexity of the problem by studying the waveforms in a 
preferred frame precessing with the binary are under way [1r39l - 
17111 . and could help in the construction of future precessing 
templates, but the detection of BHNS systems in Advanced 
LIGO remains an important challenge today (see also 072I1 for 
an updated template bank for binaries with arbitrary spins) . 
Considering the negligible influence of tidal effects on the 
waveform before the disruption of the neutron star, these is- 
sues are however better studied in the context of black-hole- 
black-hole binaries, for which longer and more accurate pre- 
cessing waveforms are available. Results obtained for these 
systems should be immediately applicable to BHNS inspirals, 
at least in the range of mass ratios considered in this work 
(tidal effects are more important for more symmetric mass ra- 
tios). 

An interesting particularity of BHNS mergers with high 
black hole spins, precession, and a high enough mass ra- 
tio that tidal disruption either does not occur or only occurs 
very close to the marginally stable orbit, is the possibility for 
the remnant black hole to receive a significant velocity kick 
from the merger. This effect is already well-known for binary 
black hole systems I73ll74ll : binary black holes with parame- 
ters similar to the BHNS mergers studied here would receive 
velocity kicks t>ki C K Si 2000 km/s (and Wkick Ss 5000 km/s 
for more favorable parameters). Kicks in BHNS mergers are 
generally much smaller, with t>kick ~ 50 km/s — 100 km/s, 
even in the precessing configurations that we previously stud- 



We continue our study of BHNS mergers at mass ratio 
q = 7, the regime currently deemed to be the most likely for 
BHNS binaries in the field. Previous results lfl8ll have shown 
that for such mass ratios, high black hole spins xbh > 0.7 
are needed for the neutron star to disrupt, even for large neu- 
tron stars (i?NS = 14.4 km). In this work, we look into 
the influence of the radius of the neutron star and the ori- 
entation of the black hole spin on the merger, focusing on 
configurations with mass ratio q = 7 and black hole spin 
Xbh = 0.9. We show that the transition between config- 
urations for which the disruption of the neutron star causes 
the formation of massive accretion disks and those where the 
neutron star just plunges into the black hole is very sensi- 
tive to the compactness of the neutron star: while a 1.4M Q 
neutron star of radius _Rns = 14.4 km leads to the for- 
mation of a massive disk (A/disk = 0.2Af Q ) and tidal tail 
(Aftaii = 0.25Af Q ), a smaller neutron star i?NS = 12.2 km 
in an otherwise identical binary forms a much less massive 
remnant (A/disk = 0.06M Q , M ta ii = O.O9Af ). For neutron 
stars with radii i?NS < 10.5 km, we expect no disruption at 
all. This indicates that in the range of stellar radii currently 
favored JUt, a black hole spin xbh ~ 0.9 is required for 
a 1.4M© neutron star to disrupt — and thus for post-merger 
electromagnetic counterparts such as SGRBs and kilonovae 
to be possible. We also note that for all but the most massive 
disks a fairly low maximum density p ~ 10 11 g/cm 3 is ob- 
served, about an order of magnitude lower than for disks of 
similar masses at lower mass ratios q = 3 — 5. This is sim- 
ply the expected geometrical effect: the radius of the disk is 



19 



roughly proportional to the mass of the final black hole. But 
that difference could significantly affect the late time evolu- 
tion of the disk: the opacity of the disk to neutrino radiation 
will be lower, and the evolution of the magnetic field is likely 
to be affected as well. 

The amount of unbound material, approximated through 
measurements of the energy of fluid elements in the limit of a 
time-independent metric, is found to be larger for these high- 
spin configurations than in previous, lower spin studies of 
BHNS mergers. The accuracy of these mass measurements 
is only ~ 50% in our general relativistic simulations — but 
for all three equations of state studied here, we find ejected 
mass M e j <; 0.01M Q . The ejecta has a velocity distribu- 
tion peaking atv/c ~ 0.2 (except for the larger neutron star, 
for which we find larger velocities which cannot be accurately 
measured at this point) and a kinetic energy Ekin <; 10 51 ergs. 
It would be a promising setup for a potential "kilonova", and 
could even be detected as a radio afterglow. It should however 
be emphasized that these massive ejecta are only produced in 
the high spin region of the parameter space, as the neutron star 
does not disrupt for low black hole spins. 

The effect of a misaligned black hole spin is also studied 
in more details. General relativistic simulations of precess- 
ing BHNS binaries had only been performed for one set of 
binary parameters before this work, for low mass, low spin 
black holes (q = 3, xbh = 0.5 l26ll ). These high spin 
configurations allow us to observe the effect of the misalign- 
ment of the black hole spin on tidal disruption more accu- 
rately. In particular, we confirm that using the radius of the 
innermost stable spherical orbit ?*isso (xbh , 0) l63ll as a way 
to predict the mass remaining outside the black hole at late 
times is reasonable: it matches the results of an aligned con- 
figuration with effective aligned spin Xefi an d nsco(Xeff) = 
nsso(XBH,0) JUt]. The qualitative features of the remnant 
can however be quite different, with the inclined configura- 
tions being slower to form a disk, and keeping more mass in 
their tidal tail. 

By extracting the gravitational-wave signal emitted by each 
of these BHNS mergers, we confirm that the effects of tides 
on the waveform of q = 7 binaries are negligible during most 
of the inspiral: tidal effects are below the numerical error in 
the simulation up to / ~ 0.5 kHz, and of the same order as the 
expected PN corrections ll3~ll HH at higher frequencies. The 
cutoff in the gravitational-wave spectrum at the frequency at 
which the neutron star disrupts is a slightly more promising 
imprint of the neutron star equation of state on the waveform. 
By combining Post-Newtonian results at low frequencies with 
our numerical waveforms, we estimate that differences in the 
neutron star radius of order 2 km could be measured in at most 
2% of the Advanced LIGO events, for a single detector at the 
current design sensitivity. Similar results can be obtained for 
the lower power lasers that are expected to be in use when 
Advanced LIGO first begins to take data, if the detector is 
tuned to observe at ~ 1.5 kHz. These estimates are however 
very optimistic, as they neglect the degeneracy between the 
neutron star radius and other binary parameters. 

Finally, we observe that some precessing BHNS binaries 
receive significant kicks from the merger, as opposed to what 



was observed for non-precessing systems. This is in agree- 
ment with results for BBH systems: the largest kicks are found 
for systems with partially misaligned black hole spins. In 
BHNS binaries, however, these results are modified by the 
fact that after the neutron star disrupts, the gravitational-wave 
signal becomes weaker. As kicks mostly arise from the emis- 
sion of gravitational waves right around merger, their magni- 
tude is significantly reduced for disrupting BHNS binaries. 
Accordingly, the largest kicks in BHNS systems are found 
for binaries with misaligned spins for which the neutron star 
does not disrupt, or disrupts very late. Thus, the kicks mea- 
sured at a mass ratio q = 7 are actually larger than for more 
symmetric systems, in opposition to BBH results. We find 
«kick ~ 345 km for our most precessing system, which rep- 
resents a lower bound on the maximum kick attainable by the 
ensemble of all similar configurations with different orbital 
phases at merger. 

One of the most important limitations of this work is that 
we do not model some critical physical effects: more realis- 
tic equations of state are required to study in detail the evo- 
lution of the tidal tail and the characteristics of the ejecta, 
and are also a prerequisite for the inclusion of neutrino ra- 
diation. Neutrinos are the main source of cooling in the disk, 
and cannot be neglected if we want to continue our evolution 
for longer than the few milliseconds after merger presented 
here. Finally, magnetic fields should also play a crucial role 
in the evolution of the accretion disk — although their evolu- 
tion requires the use of a numerical grid much finer than what 
we currently use in our simulations ll66ll . These effects are 
not expected to significantly affect the results presented here, 
which focused on the general properties of the merger and on 
the gravitational-wave signal, but they should be included in 
simulations aiming at a more detailed description of the evo- 
lution of BHNS systems after merger. 
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Appendix A: Combining PN and numerical results to assess the 
detectability of Equation of State effects in the gravitational 
waveforms 



The influence of tidal effects in the low-frequency part of 
the waveform can be estimated from simple post-Newtonian 



considerations: to leading order, the amplitude of the Fourier 
transform of the waveform is 



£>7T 2 / 3 



while its phase is 



*(/)=*o(/) + *r(/) 



(Al) 



(A2) 



where ^t(/) contains the tidal effects, and *f>o(f) all other 
contributions. In the regime considered here, 1 $>o(f) is poorly 
constrained. However, it is also identical for all 3 of our non- 
precessing configurations. Thus, we can compute the inner 
product 1 1 Sh , | for two BHNS binaries which only differ by the 
equation of state of their neutron star using hi (f) = A(f) 
and h 2 (f) = A(/)e 4 ( A *? N ^'+A0+a/) i where A*p n (/) 
is the phase difference due to tidal effects (computed from 
Eq. (T%), A(j) is an arbitrary phase shift, and a = 27r(At) 
umefty ow §dPa^9W'itrary time shift, a and A(j> are chosen to 
maximize {hi, h 2 ). Limiting our integration to frequencies 
below 0.8 kHz, we find that tidal effects during the inspiral 
are a slightly smaller effect than the disruption of the neutron 
star for the Zero Detuned noise curves (see Table II Vt. 

Considering the low-frequency and high-frequency por- 
tions of the waveforms separately result in an underestimate 
of the total \\8h\\. Indeed, each part of the waveform used dif- 
ferent time and phase shifts to maximize {hi,h 2 ). We would 



thus expect \\8h\\l ot > \\Sh\\ 2 <0 _ 8kHz 



II^H>0.8kHz- Ifwe 



assume that the tidal part of the Post-Newtonian approxima- 
tion is valid at the beginning of the simulation (which is ap- 
proximately true, as shown in Fig. [T2l . and that the simula- 
tions only differ by the effects of the neutron star equation of 
state (neglecting the residual eccentricity, as well as numer- 
ical errors), we can however obtain reasonable estimates of 
the total | \Sh\\. Indeed, we can write the waveforms from two 
numerical simulations as 

hi(f) = 4 1 (/)e ,( *» (/)+ ^ l(/)) (A3) 
h 2 (f) = A 2 (f)e i ^°W+* T ^». (A4) 

As for the Post-Newtonian expansion, the phase ^o(f) does 
not contribute to \\8h\\, while A^^ R = ^raif) ~ *t,i(/) 
can be matched to the Post-Newtonian A\&™ (/) over a given 
frequency range. Practically, we compute A'J/y R from the 
simulations, allowing for an arbitrary time and phase shift 
in one of the numerical simulations. The shifts are cho- 
sen to minimize the difference with the Post-Newtonian pre- 
dictions in the range 0.3 kHz < / < 0.8 kHz (modifying 
the matching window has a ~ 10% effect on \\Sh\\, simi- 
lar to the differences between the various PN orders). We 
then smoothly connect the amplitude and the phase of the 
Post-Newtonian and numerical waveforms through Ytot(/) = 
a(/)Y PN (/) + (1 - a(f))Y NR (f) with 



a(f) = 0.5* 1 + cos 



Af ~ fi) 

fu ft 



(A5) 



for fi < f < f u [where and /„ are the bounds of the 
matching window and Y(f) is either A(f) or A$t(/)]- We 
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also have a(f < //) = 1 and a(f > f u ) = 0. The re- 
sulting hybrids contain the information needed to estimate the 
difference | \Sh\ \ between waveforms. It is worth noting, how- 
ever, that they are not proper hybrid waveforms, and would 
be useless as templates. Indeed, we used the fact that the ony 



difference between the 3 systems considered here is the equa- 
tion of state of the neutron star to neglect the non-tidal part 
of the phase, ^o(f). However, knowledge of ^o(f) would 
be needed in order to compare an observed waveform to the 
result of our simulations. 



